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ABSTRACT 

We perform analytic linear stability analyses of an interface separating two stratified media threaded 
by a radiation flux, a configuration relevant in several astrophysical contexts. We develop a general 
framework for analyzing such systems, and obtain exact stability conditions in several limiting cases. 
In the optically thin, isothermal regime, where the discontinuity is chemical in nature (e.g. at the 
boundary of a radiation pressure-driven H il region) , radiation acts as part of an effective gravitational 
field, and instability arises if the effective gravity per unit volume toward the interface overcomes that 
away from it. In the optically thick "adiabatic" regime where the total (gas plus radiation) specific 
entropy of a Lagrangian fluid element is conserved, for example at the edge of radiation pressure-driven 
bubble around a young massive star, we show that radiation acts like a modified equation of state, 
and we derive a generalized version of the classical Rayleigh- Taylor stability condition. 



1. INTRODUCTION 

The superposition of a dense fluid above a lighter one in 
a gravitational field is p rone to the well-know n Rayleigh- 
Taylor instability (e.g. Chandrasekhar 1981): Any cor- 



rugation of the interface between them will grow expo- 
nentially, as fingers of the heavier fiuid sink in the more 
buoyant one. The Rayleigh- Taylor instability and related 
processes have found applications in various astrophys- 
ical settin gs, such as the expa nsion of supernova rem- 
nants (e.g. Ribeyre et al.||2004 ) (where inertial accelera- 
tion plays the role of the gravitational field) , the interiors 



of red giants, subject to thermohaline mixing (e.g. Char- 
bonnel & Lagarde 2010 1, or interstellar gas clouds pushed 
above the galactic plane (e.g. [Zweibel 1991) . 

One can envision several Rayleigh-'iaylor-like configu- 
rations of astrophysical interest where radiation is impor- 
tant for both energetics and dynamics. For instance, dur- 
ing massive star formation, radiation pressure overcomes 
gravity and causes the formation of bubbles of rarefied 
matter around the central star(s). Since they are overlain 
by denser infalling gas, they may be prone to Rayleigh- 
Tayl or instabilities, potenti ally aiding continued accre- 
tion (Krumholz et al. |2009 1 . Another astrophysical set- 
ting of relevance could be the interface between an H ii 
region and its neutral shell. Stellar photons are absorbed 
in the H ii region and exert a force toward the interface 
that acts like an effective gravitational field. In suffi- 
ciently dense HII regions driven by suffici ently massive 
stars this radiation force ca n be very large (Krumholz & 
Draine 2010 ), potentially destabilizing the 
Finally 
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shell ot swept-up material 



in the same vein, the 
radiation force could be significant in shaping Rayleigh- 
Taylor instabilities in supernova explosions. 

Radia tive Rayleigh- Taylor instabilitie s are not a new 
subject. Mathews & Blumenthal (1977) studied the sta- 



bility of surfaces and slabs of fully ionized plasmas and 
found instability for optically thin clouds at their far 
side and optically thick ones (using the Boussinesq ap- 
proximation) with significant am ount of neutral gas, or 
pushed at the illuminated side. Krolik (1977) studied 
the global stability of a constant-density slab under the 
Boussinesq approximation, and found, in the absence of 
gravity, instability of short-wavelength perturbations if 
radiative acceleration correlates positively with total op- 
tical depth; inclusion of gravity induced a transition back 
to the classical Rayleigh- Taylor result. 
A noteworthy related, albeit qualitatively differ ent in- 



stability was studied by Blaes & Socrates 
- hey performed a 



( 2003[ ) in the 
ocal radiative 
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optically thick regime, 
magnetohydrodynamics stability analysis of a stratified 
equilibrium, and found radiation to overstabilize acoustic 
disturbances for high enough background fiux. Radiation 
slips into rarefied regions giving rise to buoyant "photon 
bubbles". In the absence of magnetic fields, the insta- 
bility criterion requires the specific opacity to have an 
explicit dependence on the density or the temperature. 

In this study, we investigate the role of radiation in the 
linear stability of a single interface between two media, 
ignoring magnetic fields and chemical processes as well as 
the structure of the interface. We present general frame- 
works in the optically thin and optically thick regimes, 
before giving analytical solutions in limiting cases. In §2, 
we will review the fundamental equations and outline the 
model and a few generalities, while §3 applies our for- 
malism to the standard (non-radiative) Rayleigh- Taylor 
instability to illustrate how it works. In §4, we focus on 
the optically thin regime, and in particular the isother- 
mal limit, while §5 will be devoted to the optically thick 
regime, and in particular the adiabatic approximation, 
whereby the total (gas plus radiation) specific entropy 
is conserved for a Lagrangian fluid element. In §6, we 
conclude. 
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2. GENERALITIES 
2.1. Equations of Radiation Hydrodynamics 

We begin by reviewing the fundamental equations 
of radiation hydrodynamics (RHD). Beforehand, a few 
words on notation; scalars will be written in itahcs (e.g. 
a), vectors in bold (e.g. F) and higher-rank tensors in 
bold calligraphy (e.g. T). The product of two tensors U 
and V is written UV] their contraction is denoted by a 
dot for a single index {U ■ V) and a colon for two indices 
{U : V). Quantities evaluated in the frame comoving 
with the fluid will be given a subscript 0. 

In the nonrelativistic a nd inviscid limits, the RHD 
equat ions are given by (Mihalas & Weibel Mihalas 
1981: 



dp 
dt 



dt 



+ V-(pv) = (1) 
Dv 

P^ = Go-VPg+pg (2) 

-V- (?/<,v) = -PgV-v + cGO, (3) 



for the gas (mass conservation, momentum and internal 
energy equation), and: 



dEr 
1 9F 



V • F = -cG° 



G 



(4) 



(5) 



for the radiation (energy and momentum equation). Here 
p, Pg = pa^, Ug — Pg/{'-f — 1) 3,16 thc gas density, pres- 
sure and internal energy per unit volume, respectively, 
with a ~ y^ksT/m the isothermal sound speed, and 
g = — is the gravitational acceleration with (j) the 
potential. Er, F, Vr are the energy density, energy flux 
vector and pressure tensor of the radiation field. 

The rate of 4-momentum transfer from radiation to 
matter per unit space-time volume dVdt (or, minus the 
4-divergence of the radiation energy-momentum tensor) 
Q, evaluated in the comoving frame, assuming that the 
gas is in Local Thermodynamic Equilibrium (LTE), and 
that scattering is isotropic, is given by: 



Go 



Go 



/ p{KjEra - KpaT^) 



(6) 



where kj and Kp are frequency-integrated absorption 
opacity means weighted against the spectral energy dis- 
tribution (SED) of the radiation and a Planckian at the 
gas temperature T, respectively, and Kp vs. the flux mean 
(with both absorption and scattering contributions). 

When coupling of the radiation to the gas (through 
the latter term) is significant, it is useful to rewrite 
the radiation equations in terms of the comoving frame 
energy density Era — Er — 2-v ■ F/c^, radiative flux 
Fo = F — {Er + Pr)^ and radiation pressure tensor 
V r^ = Vr - (Fv vF)/c^ ( e quatio ns 95.87 and 95.88 



of Mihalas fc Weibel Mihalas (1984)): 



^+V(i?.ov-f Fo)+7'.o : Vv-H2^ 
at & 



-cGl (7) 



G - ^ ^F 



V • 7^.0 + (Fo • V) 



V • V, 



^{ErO^Vro) 



(8) 



where it is (generally) safe, for v ^ c, to drop all the 
terms containing a = Dv/Dt as well as (Fq • V) v/c^ 
and (V • v/c^)Fo. Note that this system of equations 
needs a closure, which will be obtained through various 
approximations depending on the regime considered in 
the next sections. 
If we sum equations (p|) and ([5|, we obtain the to- 



tal momentum e quatio n (equation 94.10b of Mihalas & 
Weibel Mihalas] ([l984|): 



d_ 
di 



pw 



V ■{pvv + Vr + PgIs) 



-pV0, (9) 



where X3 is the 3x3 identity matrix , and we have followed 
Mihalas & Weibel Mihalas ( 1984 ) in dropping the term 
— Ggv/c on the right-hand-side as non-dominant in flows 
with V <^ c. 

Equatio ns ([3| and ([7| can be summed to yield (equa- 
tion 16 of|Buchler|(|ig79|)): 



DEt. 



Dt 



V • Fo + Htot : Vv = 0, 



(10) 



with Etot = 

Ptot = Pg^S 



Ero + u and Htot = Etotiz + Ptot, where 
f Pr- Yet another useful form of the total 
energy equation can be obtained by adding t he scalar 
produ ct of equation ^ with v (equation 18 of Buchler 
fT979|): 



= 



d_ 
dt 



Eu 



Et, 



V • 



v + Fo + not • v^(ll 



where we have assumed the gravitational potential to be 
static. 

2.2. Model and Linear Stability Formalism 

We consider a plane-parallel background configuration, 
consisting of two semi-infinite media separated by an in- 
terface at z = 0, with medium 1 overlying medium 2 
(which one might generally think of as being more rar- 
efied). Throughout this study, we will ignore the width 
of the discontinuity, and we allow no flow accross it. The 
system is subject to a constant and uniform external 
gravitational fleld (or, equivalently, an inertial acceler- 
ation) g = —gGz and is threaded by a radiative flux F, 
which in equilibrium is independent of z and vertical. 
For the astrophysical applications considered here, both 
gravitation and radiation fields may be thought of as be- 
ing caused by a radiation source such as a massive star 
located at z — —od. 

The system of dynamical equations written in the pre- 
ceding subsection, when supplemented by equations of 
state and appropriate closures, may be cast in the form 



(12) 
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where tp is a vector of the different fields (here, phys- 
ical quantities as functions of spatial location) evolved 
in time by the (nonlinear) operator H. The equilibrium 
configuration VJcq then satisfies H{ipcq) = 0. Considering 
a perturbation Stp = ip — ^eq^ we have, to linear order 



The usual kinematic relationship 



(13) 



with dH^^^ the (linear) differential of H at -^oq- The 
problem now amounts to finding the eigenmodes of 
dH^^^, since if di?^^,((5V(0)) = 5^{t) = 

e~"^*(5'0(O). If Im(cLi) > 0, the perturbation grows and 
linear instability is declared. We therefore are interested 
in Eulcrian perturbations whose space-time dependence, 
for any quantity Q{x, z, t), is given by 



6Q{x,z,t) = 5g(z)e*('=="- 



(14) 



where the Fourier dependence in x (whereby we orient 
the axes to have k positive) is motivated by the plane 
parallel nature of the background equilibrium. Since no 
perturbed vector quantity has a component perpendicu- 
lar to both and e^, the linear problem is 2D. 

The eigenvalue problem now reduces to a set of coupled 
ODEs, supplemented by a set of relationships with no 
derivatives in z, and the former may be cast in the form: 



dz 



A{z) ■ S^, 



(15) 



where the linear operator A, in our problem, depends 
on z only through the background quantities, in turn 
completely determined by their values at z = 0* and 
the values of g and (the z component of the radiation 
flux) from the equilibrium equations. 

Since each solution to this set of ODEs corresponds to 
a set of perturbations in the space z ~ 0^, it should in 
principle be possible to analyse stability conditions as 
a function solely of quantitie s evaluated at the interface 
(rather than integrals, as in |Krolik| ( [T977| ), but he was 
considering an upper boundary tor the cloud). 

2.3. Boundary conditions 

Up to this point, nothing distinguishes our problem 
mathematically from a stability analysis of an infinite, 
single medium, be the analysis global or local in nature. 
The distinguishing characteristic of the interface prob- 
lem is the boundary conditions which select the relevant 
(z dependence of the) eigenfunctions, on which we now 
focus. 

First, we consider media that are unbounded on either 
side of the interface, so we require our modes not to blow 
up as z goes to ±oo. Thus we are focusing on "local" 
instabilities at the interface, rat her than global ones on 
a cloud scale as in Krolik (1977). This requires 



lim S'il){z) 



0. 



(16) 



We next investigate the continuity conditions at the 
interface. Let ^{x,z,t) be the Lagrangian displacement 
of the fluid element that is at position (x, z) in the unper- 
turbed state. We also denote, for any quantity Q{x,z,t) 
the Lagrangian perturbation by 



Av 



Dt 



(18) 



reduces, for the Fourier dependence adopted for our so- 
lutions and the zero-velocity background, to Sv — —iuj^. 
Since there is no fiow accross the interface, £,z{x,0,t) 
represents the vertical displacement of the boundary be- 
tween the two fluids. Thus, is continuous at the inter- 
face. 

Now consider a general flux-conservative form equation 
describing the evolution of the system: 



dm 
~dt 



Vf 



(19) 



where m is the conserved quantity, f is the correspond- 
ing flux, and s is a source term. The equations of mass 
conservation ([l]), total momentum conservation ([9]), and 
total energy conservation (111 are all manifestly of this 
form. We place ourselves m an inertial frame comov- 
ing (at time t) with the interface. (Note that for the 
general considerations we are about to make, it is imma- 
terial whether there is a net flow across it or not). Our 
purpose here is to find under which conditions the com- 
ponent of the flux /o,z (the subscript referring to the 
frame chosen) normal to the interface can be considered 
continuous accross at the interface at z = S,z. 



Integration of equation ( 19 ) accross the interface thick- 
ness yields 



fo,z{x,^z + - fo,z{x,£,z - 

5(mo) d{fo.^) 



dto 



dx 



(20) 



(17) 



where e is the thickness of the interface, and the brackets 
(...) denote averages accross the interface, i.e., for any 
function Q{x, z, t) 

(Q) = - / Q{x,z,t)dz, (21) 

which is a function of x and t. It is a consequence of the 
choice of frame and the orie nta tion of z axis normal to 
the interfac^that equation (20) has no extra "boundary 
term" . 

We expect mp, fo and s to remain bounded within the 
interface (although their z- derivatives may be large) such 
that their z-integrated averages are comparable to their 
asymptotic values on either side of the interfacej^ There- 
fore, it is already qualitatively clear that the continuity 
of fo,z will be verifled if e is "small enough" . 

^ Actually, the normal to the perturbed interface gener- 
ally differs from the z axis (defined at equilibrium), such 
that the relevant component of the flux we should consider is 

[/o,z — fo,x(d^z/dx)] /^/l + (d£,z/dx)^, but this does not differ 
from /o_z to linear order. 

^ As regards s, it is nonzero only for the momentum equation, 
where it is proportional to p and a fixed gravity; were we including 
self-gravity, we would even be able to write it as the divergence 
of a flux — (g'^/2 — gg)) /AnG so that no source term would be 
present. 
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In order to be more quantitative, we note that 

e 



(22) 



where subscripts "eq" refer to the equilibrium, unper- 
turbed value of a quantity, and for any quantity Q{x, t) 
we define 



[Q]l = Q{x,Q+,t)~Q{x,0-,t), 



(23) 



with Q{x,0^ ,t) the value taken by Q in medium 1 at 
z = (the value being extrapolated if the perturbed 
interface is actually above z = 0) and Q{x,0^,t) that 
same quantity for medium 2 (extrapolated if the per- 
turbed interface is below z — 0). We have also used 

dfo,z,eq/dz = Seq. 

Since the (scq)eoq term e ssen tially cancels (s)e in the 
right-hand-side of equation ( 20 ) , we see that the question 



of the vertical flux continuity amounts to that of A/0.2 
(which is actually w hat we will be using in the stability 
analyses). Equation (20) may be rewritten as 



dto 



dx 



(24) 



So the general condition that our perturbation must sat- 
isfy in order to have continuity of A/o,^ across the in- 
terface is that fee 5/0, X A/o,z and euiSm ^ A/o,z- If 
for example, Sfo^x Sfo^z, we obtain e ^ as might 
have been expected intuitively. 

Since application of this boundary condition to the 
mass conservation equation ([I]) does not bring any new 
information as there is no flow accross the interface, and 
since we will be making approximations to the energy 
equations, the sole important application (in this pa- 
per) of the above considerations is the z component of 
the momentum equation ([9|, where m = pv^ + F^/c^, 
f = pWzV -\- Vr ■ ^z + Pg^z, and s — ~pg. Thus we have 
/q z — V^Q + Pg, and the result will thus read 



[Afo,z]l^[SPg+6r^§-p9^.]l^0 



(25) 



(From now on, we shall drop the "eq" subscripts from 
the background quantities.) To linear order, we will al- 
ways have fo^x = 0, so the only important condition 
is eui 5m pg£,z ■ If we take Sm ^ pujS,z , one obtains 
the condition e ^ g/^'^- For lu of order the classical 
Rayleigh- Taylor result (rederived in the next section), 
this amounts to the constraint fee <C 1, i.e. that conti- 
nuity holds as long as we restrict ourselves to consid- 
ering perturbations with wavelengths much larger than 
the thickness of the interface. In the case where radia- 
tion forces are important, this is likely to be of order the 
photon mean free path or the radiation diffusion length, 
depending on the particular problem we are considering. 

3. THE CLASSICAL RAYLEIGH-TAYLOR INSTABILITY 

We illustrate the above formalism with the classical 
Rayleigh- Taylor instability, which also provides a bench- 
mark with which the upcoming results can be compared. 
In this section, we therefore ignore radiation and consider 



the two media to consist of constant-density (incom- 
pressible) fluids (as appropriate for liquids). The per- 
turbed mass conservation (here incompressibility) and 
Euler equations then read 



d_ 
dz 

-iup Sv 



Svz + ik Svx — 

d_ 

z 



dz 



-iujp 5vx + ik SPg — 0. 



(26) 

(27) 
(28) 



Solving equation (28 1 for 6 vx and recalling that Sv^ 
—iujS,z, equations (2Si and (27) yield 



with: 



dz 
A 



SP„ 



A 



SP„ 



1 (k- 



puj 







(29) 



(30) 



The matrix A here is independent of z in each medium. 
In general, for a constant 2x2 matrix A (a circumstance 
we shall encounter again) , the solution for each individual 

medium (keep in mind 5^ is not continuous accross the 
interface) may be written as 



5V(z) 



5P„ 



(31) 



where Ca,h are two constants of integration and bil)a,h are 
two linearly independent eigenvectors of the matrix A, 
with eigenvalues In this simple case, the eigenvalues 

in question are Va — k and rj, — —k. In order for ^ and 
SPg not to blow up away from the interface, it is therefore 
necessary that Ca = in the region z > and Cb = 
in the region z < 0. Sip must thus be an eigenvector of 
A, with eigenvalue —k in medium 1 (z > 0) and k in 
medium 2 (z < 0), respectively. 

To obtain the dispersion relation, we apply the bound- 
ary conditions at the interface, which in the absence of 
radiation reads 



[APg]l = [SPg ~ pgiz\\ = 0. 



(32) 



If we solve for SPg as a function of in the eigenvalue 
equation for each medium and plug into equation (32), 
we obtain 



= gk 



P2_ 
92 



Pi 
Pi 



(33) 



The instability criterion is thus pi > p2 as is well-known. 
The growth rate of the instability in the limit pi ^ p2 is 
Im(a;) = 



4. THE OPTICALLY THIN ISOTHERMAL REGIME 

4.1. Formulation of the equations 

We now consider radiation, flrst in the optically thin 
isothermal regime. By optically thin we mean that we 
can neglect attenuation and treat the radiation flux as 
constant and unperturbed in each of the two fluids, and 
by isothermal we mean that each of the fluids is kept at a 
fixed temperature via its interaction with the radiation. 
A discontinuity exists only because there is a chemical 



5 



change at the interface between the two fluids, and pos- 
sibly a frequency shift in the radiation spectrum at the 
interface as well (though the total frequency-integrated 
flux is constant). As a result, the fluid on one side of the 
interface interacts with radiation differently than fluid on 
the other side. 

One possible astrophysical realization of this situation 
is an ionization front, where fluid on one side of the in- 
terface is ionized and hot, and the radiation is domi- 
nated by ionizing photons, while fluid on the other side 
of the interface is neutral and cold, and the radiation 
there is shifted to non-ionizing frequencies. If radiation 
pressure forces dominate gas pressure ones in the ion- 
ized gas, this gas is swept into a thin atmosphere on the 
surface of the front. The downconversion of the ioniz- 
ing radiation to non-ionizing freq uencies occurs most ly 
within this t hin tra nsition region ( Krumholz & Matzner 
[2009 , prai^[20T0l ), and thus we can treat the situation 
as an interface problem 

In this limit, it is most convenient to lump the gravi- 
tational and radiation forces in equation ([2| together as: 



pg + Go = /5 g 



Hp 



(34) 



where we have ignored the difference between the co- 
moving and the reference frame, as appropriate in this 
regime in the nonrelativistic limit. If, as we shall hence- 
forth assume, the specific opacity Kp does not depend 
upon density, these two forces are exactly equivalent to 
an effective gravity field gofi = g -I- (kf/c)F = — ^offSz 
constant in each medium, but which may differ, as men- 
tioned above, between the two media. 

The equilibrium density profile on both sides of the 
interface is p cx exp (— geff^/a^), where a is the sound 
speed, so the scale height is a? /gcS- 

4.2. Stability analysis 

We now move on to the derivation of the dispersion 
relation and the instability criterion. The underlying hy- 
drodynamic equations are the same as in the classical 
Rayleigh- Taylor case, except that we replace g by gcS, 
and we relax the assumption of i ncom p ress ibility. The 
perturbed equations analogous to ( 26 1 - ( 28 1 in this case 
are 

d f d \ 

- iuj 5p + 5vz^p + p\ik6vx + ^5vz] (35) 
oz \ oz I 



-iupSvz + —SPg - gcs5p = (36) 



d_ 
dz 

—iupSvx + ikSPg — O (37) 

Eliminating Sv^ as in §3, using the isothermal equation 
of state Pg = pa^, and using the fact that dp/dz = 
— {gctf/a'^)p and for the background state, we obtain 



d 

dz 



Sp 




" iz' 

Sp 


p 




p 



(38) 



^ Strictly speaking such an interface has a flow across it, since the 
amount of ionized mass increases with time in such a configuration. 
However, for a strong D type ionization front the flux of mass 
and momentum across the ionized-neutral interface is very small 
compared to the flux reaching the front, and so we may safely 
neglect it. 



with: 



A 



\ a } 



'kay 







(39) 



The matrix A here is independent of z in each medium 
(as in §3). The general solution will thus adopt the form 
of equation (31|. Hence, we need to discuss the eigenval- 
ues of A. They satisfy the characteristic equation 



V a 



- fc^ = 



(40) 



Let us focus our attention to medium 1. For the veloc- 
ity and density perturbation to vanish for z — > -|-oo, we 
require that, for each eigenvector of A with correspond- 
ing eigenvalue Ai along which the solution has a nonzero 
projection. 



Re(Ai) < min 0, 



91 



(41) 



where gi and ai are the value of gcs and a in medium 
1. However, from the characteristic equation, we know 
that the average of the real parts of the two eigenval- 
ues is gi/2a\^ and therefore one of the eigenvalues has 
a real part that violates the above inequality, regardless 
of the sign of gi. Therefore, the solution cannot have a 
nonzero projection along this eigenvector, and must be 
an eigenvector of A. The same argument can be repeated 
in region 2 (one can e.g. change the orientation of the z 

axis to be in the exact same configuration) and thus, (5-0 
(which is not continuous at the interface) is an eigenvec- 
tor of A in each region, of eigenvalue Ai and A2. 

To derive the dispersion relation, we now intro duce 
the boundary condition at the interface. Equation (32 1 
continues to hold if we replace g with g^s, so 



a\2 



\5Pa 



pgcs 



U - 0. 



(42) 



Applying this to the characteristic equation ( 40 1 , we have 



1 



at 



a7 



.91 



1 



as 



52 



(43) 



'■2 V^S 

If we divide equation (|40 1 by A^ and equate the result- 
in medium, use of equation (43 I 



ing left-hand-sides for eac. 
yields (since Ai 7^ A2) 



A1A2 = -fc2 



(44) 



which implies the two eigenvalues have real parts of op- 
posite signs. While the constraint that the velocity per- 
turbation vanishes at -l-oo and —00 hereby reduces to 
Re(Ai) < 0, this does not guarantee that the density 
perturbation will do so in a medium where ges points 
away from the interface. In this case, one needs to fur- 
ther satisfjj^ 



(Im(..2))^ 



5cff 



Re(w2) 



> 0. 



(45) 



* We start from the inequality Re(Ai) < g\la\ (if we take 
the medium in question to be m ediu m 1), with Ai obtained from 
solving the quadratic equation ( |40| . We use the following use- 
ful relationships, holding for all complex values of z: Re(y'^) = 
^(|2| -|-Re(z))/2 and Im(v^) = sgn(Im{z))^(|2| - Re(z)) /2; the 
latter equation also defines our choice of branch cut in the complex 
plane. 
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Solv ing equation (44) for A2 and injecting into equation 
( 43 ) , one obtains tHe following quadratic equation: 



Growth rates of radiative Rayieigh-Tayior instabiiity 



^Ai 



a, = 0, 



(46) 



with h — 5ia2 — .92 'ii: from which one deduces that 
sgn(Re(Ai)) = sgn ((pigi - p2ff2)Re(a;2)) ^^j^j^ ^^^^ ^ 

p(0*)). Therefore, if pigi > P252, 1^ cannot be real and 
the configuration is unstable. Indeed, even if Im(cj) < 0, 
we sir nply nee d to take the complex conjugate of equa- 
tions (38)-(39): we then see that Si/;* corresponds to a 
perturbation (still satisfying the boundary conditions) 
with the same wavenumber k but with complex frequency 
uj* , and which is therefore a growing eigenmode. In a 
more coordinate-free manner (since upper and lower are 
not well-defined if the direction of the effective gravity 
switches sign across the interface), the instability arises 
if the effective weight per unit volume toward the inter- 
face overcomes that away from it. A corollary is that, 
should goff point toward the interface in both regions, 
the equilibrium is unequivocally unstable. 

In order to completely prove the sufficiency of the cri- 
terion, we need to show that Re(a;^) < (which ensures 
the inequality (45)) is actually allowed by the dispersion 
relation. To do so, we first solve for Ai from equation 



(46) after eliminating the quadratic term with equation 
(5D). We obtain 



Ai 



(fca) 



hk^' 



(47) 



with a = \J a\ +02- Since the same formula holds for A2 
if the subscripts "1" and "2" exchange roles (and thus h 
switches sign), the equation (44) yields the desired dis- 
persion relation 

= 0;* - Kkafuj"^ -t- {(kaf -t- 5132^^] w'' 

+ k\g^-g^)hu? -k'^h', (48) 

which is fourth order in cj^. As this polynomial al- 
ways has one negative real root (in terms of w^), pro- 
vided kh ^ Q, the sufficiency of the instability criterion 
piQi > /92<?2 (or, equivalently, /i > 0) is proven. The 
polynomial has two real roots, only one of which is phys- 
ically allowed (the other being opposite to leading order) , 
asymptotically given by 



2 _ hk h{gi+g2){al - aj) . ^ / 1 

^ — t" I 7r~R ' \ T 

^2 2a^ \ k 



, Pi.9i - P2.92 

Pl+ P2 

(.91 + 32) {P2 - Pi) (Pi.9i - P252) 



2a2 (P1+P2) 



+ 



(49) 



In the long- wavelength limit, if the instability crite- 
rion is satisfied, oj^ is given by (/ifc)2min (1/gi, —1/(72) 
if gig2 > and —y/—gig2k if 51(72 < 0. The instability is 
a "pure" instability (in the sense that oj is purely imag- 
inary). Figure [1] shows a calculation of the growth rate 




Fig. 1. — Growth rate s = Im(tj) of the isothermal, optically-thin 
radiative Rayleigh-Taylorinstability, calculated numerically from 
the dispersion relation ( |48| . The value shown corresponds to the 
fastest growing mode, i.e. to the largest value of s. Growth rates 
and wavenumbers are nondimensionalized through combinations 
of a = (a^ -I- al)-*^/^ and g = (g^ + 02)^^^ t E'nd as such depend on 
two dimensionless parameters I^g/g = {gi — g2)/g and the Atwood 
number yl = {pi—p2)/{p\+P2) (plus sign information, sgn(gi+g2)) 
which are bounded by 2^/^ and 1 in absolute value, respectively. 
Here, we have fixed A = 0.5 (and gi + g2 > 0) and varied l\g/g 
with values (from top to bottom) 2^/^, 1, 0.5, and —0.5 



as a function of various parameters. Note that equation 
( 49 ) agrees with earlier results for compressible Ray leigh 



ay lor instab ility without radiation (equation 23 of Shiv- 
amoggi||2008 1 if we take gi = 32 ■ 



4.3. Sample application: radiation pressure- driven HII 

regions 

We conclude this section with a sample application for 
the case of the ionization front around an H 11 region 
where radiation pressure signifi c;antly affects the d ynam- 
ics (e.g. the 30 Doradus region, Lopez et al.||2010 ). Con- 
sider such a region powered by a star cluster ot luminosity 
expanding into a uniform ambient medium of num- 
ber density n, sweeping up ambient gas as it expands. 
We will neglect the gravitational pull of the star cluster, 
which is significant only early in the evolution. During 
the radiation-dominated phase of the expansion, which 
applies when the H i l radius r <g; tq, the radius of the 
H II after a time t is ( Krumholz k Matzner||2009 ) 



ra{t/to 



,1/2 



(50) 



where tq = lliy pc, to = hQ^' nj Myr, = L^/IO"^ 
Lq, ng = rt/10^ H nuclei cm"'^, and we for all other quan- 
tities have ad opted the fiducial parameters of |Krumholz| 
Ifc MatznejPl for the embedded case. Our choice of lumi- 
nosity and density are motived by the example of the 30 
Doradus H 11 region, which is driven by a central cluster 
of luminosity 1.7 x 10^ Lq. 

The acceleration of the shell, and thus the effective 
gravitational force toward the front in the frame comov- 

^ except that we take 1/1 = 3.3, and we corre( :t a factor of 2.2 
error in equation 4 of Krumholz & Matzner - see [Fall et al.| | |2010[ l 
for details 
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ing with the front, is 



9i 



( t 



-3/2 



-1.1 X lO^^^if^ng"^ 



-3/2 



cm s ^, (51) 



where positive sign corresponds to goff' pointing toward 
the cluster. Note that, for simpUcity we have assumed 
that the material outside the transition region near the 
front is optically thin to the (non-ionizing) radiation that 
emerges from inside it. Within the shell, the opacity is 
dominated by dust (except near the transition region, 
where the neutral fraction is high enough for neutrals to 
contribute significantly). Radiation exerts a force toward 
the front, inducing an outward force per unit mass in the 
frame of the front given by 



92= gi 



KpL 



= -8.8 X lO^^Kgif^ — 



-1/2 



cm s 



(52) 



(53) 



where K3 — cm^ g^^j ^nd in the numerical eval- 

uation we have dropped gi since it is small compared to 
g2- The normalization of Kp is chosen because 10'^ 
cm^ g~^ is a typical dust opa city for radia tion at the 
color temperature of an O star ( Draine|[2003 ). 



If we adopt sound speeds of ai = 0.19 km s~ and 
a2 = 9.2 km in the neutral and ionized gas (appro- 
priate for molecular gas at 10 K and ionized gas at 7000 
K, respectively), pigi - P252 = (pi - P2)5i + P2HpF/c w 
P2KpF/c > 0, and we find that this configuration is un- 
stable. The inertial force term pigi points away from the 
interface and is therefore stabilizing, and instability oc- 
curs only because it is overcome by the larger ^232 term 
that is dominated by radiation force. In the short wave- 
length limit, which applies for all perturbations sm aller 
than ro, the growth rate is given by (using equation 49) 

^/hk ai J — - 

im(a;) w w — \J gik. 

a 02 

Plugging in our fiducial values, we have 



(54) 



Im(a;) « 0.08 



^ 1/2 / , \ 1/4 , , ,0 
K3 \ lto\ /ro\ 1/2 



X 



Myi-\ (55) 



and we learn that modes with A/ro < 0.01 (A < 0.1 pc 
for our fiducial parameters) will be able to grown signif- 
icantly in the few Myr lifetime of the stars driving the 
H II region. This may explain the small-scale filamen- 
tary structures seen around the edges of 30 Doradus and 
similar radiatively-driven H 11 regions. 

5. THE "ADIABATIC" RAYLEIGH-TAYLOR INSTABILITY 

5.1. Formulation of the equations 

We now consider the stability of an interface where 
both sides are optically thick (although we discuss how to 
relax this requirement for the lower medium below) and 
in radiative equilibrium. In this regime, the radiation 
field is a Planckian locked at the gas temperature and 
the comoving frame pressure tensor may be taken to be 



isotropic (scalar) and given by Vro = (i?ro/3)l3, with 
Ero = Equation (|8| may then be approximated 

by 

Go « -yPro. (56) 

which can be lumped with the gas pressure force. In 
making this assumption, we require that the photon 
mean free path l/{Kpp) be smaller than the wavelengths 
of the perturbation and the characteristic lengthscale 
of variation of the background equilibrium. The lat- 
ter may be defined as L = imn (Ptot / pgi Proc/ KpF) = 
min (1 -|- X, x/E) a? j g, where for convenience we define 



El 



KpFo 

gc 

PrO 



(57) 
(58) 



Physically, E measures the Eddington ratio of the back- 
ground state (i.e. the ratio of radiation force to gravita- 
tional force), while x measures the relative importance 
of radiation and gas pressure. 

Considerable simplification of the problem is achieved 
if we are allowed to drop V • (5Fo in the energy equation. 
(Appendix [A| shows the system of equations without this 
approximation.) Indeed, the energy equation (10 1 can be 
rewritten as 



Ds 
^ Dt 



(59) 



with s the specific entropy of the gas plus radiation fiuid. 



m(7 - 1) ^'^ 



PT 



Perturbation of equation ( 59 ) yields 



~Dt 



V • (5Fo 

pT 



(60) 



(61) 



If we can actually disregard the right-hand-side, we ob- 
tain 

As = 0. (62) 

Such an approximation, henceforth referred to as the 
"adiabatic approximation" , holds in the limit of high 
optical thickness. More precisely, when discussing the 
validity of the upcoming calculation, we will require 



5Fo < EtotSv. 



(63) 



A possible astrophysical realization of this configuration 
is the wall of a ra refied bubble blown by a massive star 
in formation (e.g. Krumholz et al.|[2009 ). The system is 
close to the Eddington limit, so the gravity and lumi- 
nosity of the central star (at z — —00) nearly balance. 
Medium 1 would correspond to dust-laden gas infalling 
from the protostellar core while medium 2 would refer to 
the rarefied bubble. Flow of matter across the interface 
is very slow and can therefore be neglected. 

5.2. Stability analysis 
While the perturbed mass conservation equation is the 



same as in §4 (equation (35 1), the two components of the 
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perturbed Euler equation read 

d 



iujpSv;, ^-^SPtot - gSp (64) 
oz 

-iupSv^ = -ikSPtot (65) 
The total pressure continuity at the interface is 



(66) 



We ehminate Sv^ through equation (651 as previously 
and Sp through the following formula (from the expres- 
sion of s): 



Sp 
P 



c- 



.6Pt. 



Pn 



1 + 4a; TO 
D k, 



-6s 



(67) 



combined with equation (62). Here 



— [ 12x+ 
D\ 7-1 



D^Wx^ + 20x - 



1 



7-1 



(69) 



The ODE system in z then reduces to a 2 x 2 matrix 
A, defined by 



d_ 
dz 



6 



= A 



6 
5Ptot 



with 



A 



where 



(70) 



(71) 



B= — [ 1&{E - l)x'^ + (24:E - 8)x 



E{5 



7-1' 



1 



4(7 - l)x 



(72) 



Let us restrict attention to cigenmodes actually lo- 
calized at the interface, i.e., which vanish on a vertical 
lengthscale small compared to L. This we will refer to as 
the "evanescence condition" . Then, the matrix A above 
may be treated as a constant in each region (we pro- 
visionally consider the adiabatic approximation to hold 
in the lower region too). Since it is traceless (and thus 

the two eigenvalues are equal and opposite), -0 niust be 
an eigenvector of A (as it was in the previous section) 
in each region, of eigenvalue Ai and A2 in media 1 and 
2, where the eigenvalue have negative and positive real 
parts, respectively. 

Combining the first row of the eigenvalue equation for 
— 0^ and the continuity of the pressure, one finds that 



\{uj/k) — g 
'1 - C(w/fca)^ 



0. 



(73) 



More explicitly, the dispersion relation implied by the 



above reads 



Pi 



1 - Ci {uj/kaiY 



Ci 



Ci(Bi -Ci).g2 



fc2 + Si 



kg 



P2 



C2 



C2iB2-C2)g^ 



kg 
a2Uj 



(74) 



This equation may (in principle) be manipulated to 
achieve polynomial form, but with a degree of 20, and 
a loss of sign information. 

We now further specialize to the case where p2 <C pi- 
Equation ( 73 ) then reduces to 



(75) 



This equation we would have obtained more directly had 
we set AiPtot = A2Ptot — and thus its validity is 
not endangered if the adiabatic approximation is violated 
in medium 2, or if the latter is optically thin, provided 
we consider this rarefied medium to have an imposed 
radiation field. We observe that APtot = combined 
with As = entails AT = and Ap = at the interface. 



From equation (75), one deduces that we must have 
Re(aj2) < if g > 0, in order for the perturbation not 
to blow up for z — > -l-oo, and thus, we have instability. 
Dropping henceforth the '1' subscripts, the evanescence 
condition then translates to 



k 



< a\ min ( 1 



E 



(76) 



The dispersion relatio n fo llows from setting the right- 
hand-side of equation (74) to zero: 



C 



, , fC{C-B)g^ 



Vfc' = 0, 
(77) 

which is third order in uP' and always has a negative root 
(in terms of w^), whence the instability qualifies as a pure 
instability. 

In the long- wavelength limit, the complex frequency 
converges toward a finite value given by 




Ax 



4(7-l)a; 



16a;2 + 20a; 



7 
7-1 



(78) 



which amounts to a Brunt- Vai'sala frequency modified 
by compressibility, corresponding to a positive growth at 
(and even somewhat below) the Eddington limit. Indeed, 
the background entropy gradient 



ds 
dz 



T 



l + Ax~E 



Ax 



7 

4(7 - l)x 



(79) 
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is already negative at the Eddington limit, since a large 
radiative flux corresponds to a large temperature gradi- 
ent in the optically thick limilj^ However, as we shall 
see in the next section, this long-wavelength limit often 
violates the evanescence condition (except for a nonzero 
range if B min (1 -I- x, x/ E) ^ 1), in which case this con- 
vective insta bili ty, although physically sensible because 
of equation (791, is not captured quantitatively by the 
present calculation. 

Since the evanescence condition allows one to neglect 
the highest order term and the first ter m in side the 
parentheses of the second one in equation ( |77[ ), we can 
provide the following more explicit (approximate) for- 
mula: 



Growth rates of radiative Rayieigh-Tayior instabiiity 




(80) 



For large wavenumbers, this coincides with the standard 
Rayleigh- Taylor result, and indeed, physically, the in- 
stability in question is really Rayleigh- Taylor in essence: 
this results from the fact that, in the adiabatic ap- 
proximation, the system is completely mathematically 
equivalent to a single fluid with a modified equation 
of state (total pressure law). We emphasize that this 
"adiabatic Rayleigh- Taylor instability", as we will call 
it for ease of reference, is not an interface counter- 
part of the loc al RHD instability studied by |Blaes fc 
Socrates (20031, where acoustic disturbances are over- 
stabie, whereas perturbations are here incompressible at 
the interface. Mo reover, the instability criterion of Blaes 
& Socrates (20031 involves the opacity law, while opacity 
is actually absent in the equations under the adiabatic 
approximation, except as a constraint for their validity. 
The domain of validity of the above calculation we shall 
now discuss. 

5.3. Validity of the adiabatic approximation 
_Since SPro = —^zidPro/dz) — Epg^z, the requirement 



(63), which amounts to kcSPro/np <C EtotSvz, is equiv- 
alent to 

k Etot 

Figure [2] shows the growth rate calculated for example 
parameters, with regions forbidden by the evanescence 
and the adiabatic constraints shaded. As mentioned 
previously, we see that the long-wavelength limit vio- 
lates the evanescence condition, while at large enough 
wavenumbers {k > g (£'tot/-Fo)^ min (1 + x/i?) if we 
take « ~9k), the adiabatic approximation breaks 
down. 

Constraints (|76|) and (81) are compatible (i.e., allow a 



range of wavenumbers in which the above calculation is 
valid) if 



Fo < i^totaA/min 1 



E 



(82) 



Inasmuch as it is a diffusion term that becomes impor- 
tant (in the energy equation) as the adiabatic approx- 
imation breaks down, it may be suspected that larger 

Since Fq < Eroc, we must have E/x < KpPg/g, and by a large 
margin at that if the diffusion approximation is to hold. 



Not local ^^.^fi^'^'^'^ 


E=100 x=l&0 




E=ix=i ^„,:f<-y^ 






/ 




/ 




Not adiabatic 



Fig. 2. — Growth rate s = Im(aj) of adiabatic Rayleigh- Taylor 
instability, calculated numerically from the full dispersion relation 
njj (solid line) with the classic Rayleigh- Taylor result overplotted 
(dotted line). Shaded in yellow is the region of the k-s plane where 
the evanescence constraint (the decay of perturbations on length- 
scales smaller than the background equilibrium scale height) is vi- 
olated, and in red that of breakdown of the adiabatic approxima- 
tion (conservation of specific entropy for gas plus radiation) breaks 
down, for a given value of Fo/Etoti (which is here arbitrary, since 
the calculation does not depend on _Fo). For completeness, we note 
that allowable growth rates and wavenumbers are also bounded by 
radiative equilibrium and optical thickness requirements, respec- 
tively. We have plotted the growth rates for three representative 
sets of parameters E = 0, x = 1 (for which the growth rate vanish 
in the long wavelength limit) E = x = 1 and E = x = 100 (for 
which the growth rate converges toward a finite value in the long 
wavelength limit). In the short wavelength limit, the growth rates 
are identical to the classical Rayleigh- Taylor value. 

wavenumber disturbances are damped (yielding a max- 
imum growth of order -^/min (1 + x, x/ E)gE^^at/ Fq — 
^min(l + x,x/E)nEtot / Ec) . 

5.4. Sample application: bubbles around massive 
protostars 

As a sample application, we consider radiation-driven 
bubbles aroun d massive stars, such as tho se seen in the 



simula tions of Krumholz et al. (2009 



(2010). In these simulations, the radiation emitted by 



or |Kuiper et al 



an accreting massive star exerts a force stronger than 
gravity on the dusty gas around it that is trying to ac- 
crete. Above and below the midplane of the dense ac- 
cretion disk, the radiation drives an expanding shell of 
material into the surrounding protostellar core. Inside 
the shell is an evacuated low-density region filled by ra- 
diation, while outside there is dense swept-up dusty gas. 
We are interested in exploring the stability of the inter- 
face between the low-density bubble and the high density 
shell, where the flow is near the Eddington limit. Given 
the relativel y sharp edges of the radiation bubbles that 
form in the iKrumholz et al. (2009) simulations, it does 



appear that they lend themselves to an interface stabil- 
ity analysis such as the one we present here. The con- 
tinuous medium case treated by Blaes & Socrates ( 2003 1 
ind eed does not see m to be of relevance here, because 
the IKrumholz et al.1 simulations — which did not include 
magnetic fields — do not satisfy a priori their local hy- 
drodynamical instability criterion (equation 58, or 63 if 
one takes gas and radiation temperature to be equal) 
as their implemented specific Rosseland mean opacity is 
independent of gas density. 
Consider a star of mass A/* and luminosity L* that 
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log A [AU] 



Fig. 3. — Radiation Rayleigh- Taylor instability growth rate for 
the bubble around a massive star, as a function of the wavelength 
A = 27r/fc of the perturbation, computed for the parameters Ts = 1, 
L5 = 1, p_i6 = 1, r4 = 1, and M = 100 Mq. The blue line shows 
the computed growth rate, while the pink and yellow regions indi- 
cate where the conditions for the in stab ility to be adiabatic (equa- 
tion^!} and evanescence (equation |76|l conditions are violated. 




Fig. 4. — Same as Figure I 
is = 0.1, P— 16 = 1, r4 = 1, ari3 



log A [AU] 



but for the parameters = 1, 
M = 10 Mq . 



inflates a bubble of radiation to a steflocentric distance 
R. The material in the bubble wall has a density p and 
a temperature T. The mean mass per particle is m = 
2.33 mjj and the gas has 7 = 7/5, appropriate for warm 
molecular hydrogen. The bubble wall is defined by the 
condition E \, and we also have 



(83) 



where Ts = T/llOO K (where 1100 K is the dust sublima- 
tion temperature in the simulations, and thus the tem- 
perature at the edge of the bubble) and p-ie = p/10~^^ g 
cm~^, a typical density in the bubble wall in the simula- 
tions. Fir st w e can check that the compatibility condition 
(equation 82 1 for application of the adiabatic approxima- 
tion is satisfied. Doing so, we find 



50 



T 1/2 ' 



(84) 



where r4 = r/10'' AU, L5 = L/10^ Lq, and we have cho- 
sen normalizations for our parameters based on typical 
bubble properties seen in the simulations. Since adia- 
baticity requires that this ratio be 3> 1, we see that the 
condition is satisfied, and instability is guaranteed since 
P2 < Pi 



To obtai n th e growth rate, we plug into equations ( 68 ) 
(691), and ^ for B, C, and D, and using g = GM.JW 



numerically solve the dispersion relation ( 77 1 to find 



negative value of w^). We do so for two example sets 
of parameters in Figures |3] and [4j In each case we see 
that the instability growth time is below 1 kyr, short 
compared to the 100 kyr formation timescale for the 
star. At small wavelengths the instability is likely to be 
suppressed by diffusion, sin ce the solution violates the 



adiabaticity constraint (81 1, but at large wavelengths 
(in some cases comparable to the physical size of the 
bubble) the constraint is satisfied and instability occurs. 
This instability explains the behavior observed in the 
[Krumholz et al. [ simulations, where large modes grow un- 
stable and allow accretion onto massive stars that are for- 
mally super-Eddington. Although we have argued above 
that photon bubble instabilities were not relevant to in- 
terpret these simulations, this should not necessarily be 
extended to the "real-world" problem of massive star for- 
mation, e.g. because we have ignored magnetic fields 
(which have been detected observationally) , un like the 
local s i mulat ions of circumstellar enveloppes by |Turner| 



et al. 



( |2007[ ). The relative importance of the photon 
bubble and the radiative Rayleigh- T aylor instabilities (as 
well as other processes, see e.g. |Zinnecker fc Yorke|2007 1 



for massive star formation has yet to be determined. 
6. SUMMARY AND CONCLUSION 

We have studied the linear stability of a plane parallel 
superposition of two media separated by a thin interface, 
with both gravity and radiation force, and given results 
for two analytically tractable limiting cases. In these two 
cases, the role of radiation in these Rayleigh- Taylor-like 
instabilities is qualitatively different. 

In the optically thin, isothermal limit, with a constant 
flux and a constant specific opacity in each medium, as- 
sumed to be chemically distinct, radiation acts like an 
effective gravitational field, which generally is different 
on either side of the interface. Linear instability occurs 
if the effective gravity per unit volume toward the in- 
terface overcomes that away from it, which in the case 
of a continuous effective gravity reduces to the ordinary 
Rayleigh- Taylor criterion on the Atwood number. This 
instability might contribute to the asymmetry of H 11 
regions. 

In the opposite limit, if the upper medium is optically 
thick and satisfies the approximation that the total spe- 
cific entropy of the gas plus radiation fluid is conserved, 
assuming the lower medium to be rarefied, one finds that 
perturbations that vanish away from the interface more 
rapidly than the background equilibrium scale height are 
unstable. In the short-wavelength limit, the instability 
is indistinguishable from the classical Rayleigh-Taylor re- 
sult, since the adiabatic approximation reduces the sys- 
tem to a single fluid, where the radiation force is part 
of the pressure force. Sufficiently close to the Eddington 
limit, the growth rate converges toward a finite value 
in the long-wavelength limit because of the negative en- 
tropy gradient. This instability could pertain to massive 
star formation by accre tion beyond the Eddington limit 
(|Krumholz et alj[2009|). 



the fastest growing mode (i.e. the one with the largest 



Yet other regimes not studied in this work could be 
qualitatively different. For example, it is conceivable 
that the local radiative hydrod ynamic overstabilities in 
optically thick media studied by Blaes & Socrates ( 2003 1 
have an interface counterpart, where radiation slips into 
underdense regions near the interface. Those instabil- 
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ities are not captured by the adiabatic approximation, 
and investigation of how the instabihty behaves beyond 
the adiabatic hmit is left for future work. 
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APPENDIX 

THE MATRIX A IN THE OPTICALLY THICK REGIME WITHOUT THE ADIABATIC APPROXIMATION 



The relevant perturbation equations are those mentioned in §4.2.1, plus the perturbed energy equation, diffusion 
approximation closure, and radiative flux continuity: 



d 
oz 



5¥o 



X 

with X = Kpp- The matrix A is 4 x 4 and is deflned by: 



V • (5Fo + fftotV ■5v = 
'-VSPro - Fo51nx 



0, 



d 
dz 







L' 


SPtot 




SPtot 


SPrO 




SPrO 






. SFo, _ 



(Al) 
(A2) 
(A3) 

(A4) 



and given by: 







Pa 



C Pg 

Pa 



withe„=§p^, andeT=§^, . 

P amp y ^ olnT \p 



Pa ^ iPrO 



(l + 0p)(i^ 



APrO 

4- 

^ 4P, 
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